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ABSTRACT 

Master's  Thesis  which  discusses  nature  of  allocation 
problems.   Danskin  Algorithm  for  solution  of  a  convex  func- 
tion to  be  minimized  over  a  closed  convex  set  is  developed 
An  example  of  application  involving  solution  of  a  3600 
variable  allocation  problem  using  a  computer  is  provided. 
Paper  includes  analysis  of  solution  and  discussion  of  pro- 
blems encountered  in  computer  application.. 
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I.   INTRODUCTION 

A  large  selection  of  problems  faced  within  the  mili- 
tary and  commercial  environments  consist  of  finding  the 
"best"  way  of  using  available  resources.   In  most  cases 
the  "best"  method  of  using  these  resources  is  that  which 
provides  the  most  desirous  return  to  the  using  agency  or 
organization.   Typically,  this  return  can  be  represented 
as  a  payoff  associated  with  the  use  of  resources,,  or  in 
mathematical  terms  it  is  a  function  -  the  value  of  which 
is  dependent  on  and  determined  by  the  amount  or  value  of 
resources  expended.   Naturally,  the  availability  and  use 
of  these  resources  are  bounded  and,  therefore,  act  as  a 
constraint  on  the  using  organization.   This  constraint 
could  be  that  the  organization  is  required  to  use  some  set 
minimum  amount  of  a  particular  resource,  or  in  the  more 
usual  sense  only  a  given  maximum  amount  of  the  resource 
is  available  for  use.   In  the  former  case  the  resource  is 
termed  lower  bounded,  and  in  the  latter  case  the  resource 
is  considered  upper  bounded.   In  many  such  problems  the 
resource  is  both  upper  and  lower  bounded  and  usually  some 
allocation  of  the  resource  between  the  bounds  is  acceptable. 

If  the  situation  exists  where  the  bounds  on  a  resource 
are  known  but  where  varying  amounts  of  the  resource  can  be 
expended  in  a  number  of  different  ways,  then  the  manner  in 
which  the  resource  is  allocated  becomes  of  critical  importance 
to  the  allocating  organization  in  terms  of  the  payoff  it 


receives  from  choosing  a  particular  allocation  scheme. 
By  simple  use  of  a  crystal  ball,  random  selection  device, 
or  some  other  means;  and  by  exercising  care  not  to  violate 
the  bounds,  the  organization  can  adequately  allocate.   How- 
ever, it  cannot  be  sure  it  has  found  the  best  allocation  » 
scheme . 

If  the  payoff  realized  can  be  expressed  as  a  mathe- 
matical function  of  the  resources  allocated,,  then  the  pro- 
blem becomes  susceptible  to  solution  by  mathematical 
programming  methods.   In  many  cases  then  the  "best"  or 
optimal  allocation  of  resources  can  be  found  and  proven  to 
be  optimal.   These  type  problems  can  be  described  as  "al- 
location" type  problems.   Typically,  the  value  representa- 
tion of  a  resource  variable  is  a  non-negative  number  less 
than  one  and  indicating  the  percent  of  the  total  resource 
available  which  is  to  be  expended  in  the  particular  manner 
or  activity  associated  with  the  given  variable. 

The  formulation  of  the  payoff  function  is  usually  the 
deciding  factor  in  determining  what  type  of  mathematical 
programming  method  could  be  used  to  find  the  optimal  al- 
location or  optimal  solution  to  the  problem.   For  instance, 
if  the  payoff  function  is  linear  in  nature,  then  it  may  be 
possible  to  solve  the  problem  by  using  existing  linear  pro- 
gramming methods.   However,  many  payoff  functions  cannot  be 
expressed  as  linear  relationships  even  when  attempts  are 
made  at  reasonable  approximations.   In  this  case  one  must 
turn  to  other  methods.   In  some  cases  it  may  be  true  that 


no  known  mathematical  method  for  optimizing  exists.  In  this 
case  it  may  be  proper  to  reformulate  the  problem  or  possibly 
return  to  the  crystal  ball. 

This  paper  will  address  a  particular  category  of  al- 
location problems,  develop  a  method  for  solving  problems 
which  fit  into  this  category,  and  provide  a  detailed  exam- 
ple of  a  practical  problem  solved  by  the  method  developed. 
Credit  for  development  of  this  method  belongs  to  Dr.  John 
M.  Danskin,  Professor  -  U.  S.  Naval  Postgraduate  School. 
This  method,  employing  an  algorithm  termed  the  Direction 
Finding  Algorithm,  was  presented  by  Dr.  Danskin  in  various 
classroom  lectures  presented  during  his  tenure  at  that 
school . 

A.    PRESENTATION  OF  THE  PROBLEM 

Consider  an  allocation  type  problem  in  which  the  pay- 
off function  is  convex  (or  concave)  and  it  is  decided  that 
the  optimal  solution  is  that  solution  which  minimizes  (or 
maximizes)  the  payoff  function.   Now,  consider  possible  al- 
location of  resources  to  be  represented  as  X  =  (xi,x2,...x  ) 
where  xx  =  percent  of  resources  expended  in  activity  #1, 
X2  =  percent  of  resources  expended  in  activity  #2,  etc.;  for 
i  =  1 ,  .  .  .  ,n . 

As  indicated  earlier,  certain  constraints  must  be  heeded 
in  finding  the  optimal  solution.   These  constraints  can  be 
formulated  as  follows: 

(1)   It  is  necessary  and  desirable  to  use  up  all  of  the 
resources  available, 


(2)  It  may  or  may  not  be  necessary  to  use  some  per- 
cent of  the  resource  in  any  given  activity, 

(3)  The  percent  of  resource  expended  in  activity  i 
cannot  exceed  its  upper  or  lower  bounds,  and 

(4)  Negative  use  of  the  resource  is  meaningless  and 
not  possible. 

This  problem  as  described  can  then  be  mathematically 
expressed  as : 

maximize  f(x)  where  f(x)  is  concave 
or 

minimize  f(x)  where  f(x)  is  convex 

subject  to  constraints: 

E  x.  =  1 

i 
l 

a.  <  x.  <  b. 
l  -   l  -   i 

where  a.  and  b.  represent  the  upper  and  lower  bounds  re- 
spectively of  the  allocation  to  the  itn  activity.   It  then 
follows  that  the  following  relationships  are  true: 

0  <  a.  <  b. 
i    i 

Z  a.  <  1  <  E  b.  . 
i   x       i   * 

Essential  to  the  development  and  use  of  the  Direction 
Finding  Algorithm  is  the  fact  that  the  payoff  or  objective 
function  must  be  dif ferentiable .   Therefore,  this  require- 
ment that  the  gradient  exists  and  can  be  calculated  is  con- 
sidered as  an  added  constraint. 


This  category  of  allocation  problems  just  described  is 
the  category  for  which  the  about- to-be-developed  Direction 
Finding  Algorithm  can  be  used  to  find  the  optimal  solution. 


II.   DERIVATION  AND  DEVELOPMENT  OF  THE  ALGORITHM 

In  this  section  the  Direction  Finding  Algorithm  will 
be  developed  and  a  step-by-step  procedure  described  for 
its  application. 

Recall  the  original  problem  is  to  maximize  (minimize) 

a  concave  (convex)  dif ferentiable  function  of  the  vector 

X  =  (x1,x2,...,x  )  with  constraints  as  follows: 

Z   x.  =  1  0  <  a.  <  b- 

1  -   l     l 

l 

a.<x.<b.         £   a.  <  1  <  £  b.  . 

i-i-i        .    i        .   i 

The  basis  of  the  algorithm  is  to  find  the  direction 
of  fastest  increase.   This  means  it  is  necessary  to  maxi- 
mize the  Directional  Derivative  which  is  defined  in  Ref.  [1] 
as 

D   f(x)  =  E  v.  f   (x)  =  I  fi-Y- 

Y    K     J  .   '  1    X.  K     J  .         l'l 

'        111 

where  Q    represents  the  gradient  and  y  represents  the  direc- 
tion of  the  gradient. 

To  develop  the  algorithm,  the  following  constraints  are 
imposed  on  y,  the  direction  vector: 

£  Y-  =  0 
-  '  l 

l 

E  7?  -  1 
i   x 

Y •  >  0  if  x.  =  a. 
l  -       l     l 

Y-  <  0  if  x.  =  b. 
'  l  -       l     l 


The  entire  problem  then  becomes 

maximize   Eft*y 

subject  to   E  y2    =  1 

E  y.  =  0 
1 

Y-   >  0  if  x-  =  a. 
1 1  l     l 

Y-   <  0  if  x.  =  b. 
1  i   -       i     i 

with  the  assumption  that  the  maximum  is  non-negative. 

The  first  step  to  solving  this  problem  will  be  to  show 

necessary  conditions  for  solution  by  application  of  the 

well-known  Kuhn-Tucker  (K-T)  Conditions  which  can  be  found 

in  Ref.  [2]  and  many  other  places  in  recent  literature  on 

non-linear  programming.   The  constraints  can  be  rewritten 

and  LaGrange  Multipliers  assigned  as  follows: 

CONSTRAINT  MULTIPLIER 

E  Y?  -  1  <  0 

'  l 

1  -  E  y2.    <  0 

I  i  - 

E  Yi  <    0 
-E  Y;L  <    0 

- Y .  <  0  if  x.  -  a. 

I I  -       i     l 

y.  <  0  if  x.  =  b. 

'  i  -       i    i 

Application  of  the  Kuhn-Tucker  Theorem  then  results  in 
V  (ft'Y)  =  X'V(E  Y?"l)  +  A"V(1-Zy?) 
+  U'V(E  Yi)   +  u"V(-E  Yi) 
+  E'C-v^VCy^  +  E-'Ctt^VCyP 
10 


X! 

l 

i 

"i 

u|' 

1 

TT  . 

1 

where 


V   as  usual  represents  the  gradient 

E'   implies  summation  over  those  i's  for  which 


x.  =  a. 
1    1 

E'1  implies  summation  over  those  i's  for  which 

x.  =  b - . 
1    l 

If  at  this  point  only  the  ith  component  of  the  gradient 
is  considered,  the  above  equation  leads  to 

n .  =  y  ?  ( 2  r  •'  +  2  y !  ' )  +  y  •  -  u ! '  -  v ..  +  it. 
1    '  1 v  '  1     1   J  1    1     1    1 

where  again  v.  enters  consideration  only  if  x.  =  a.  and  ir. 
6     i  7      l     l      l 

is  considered  only  if  x.  =  b-.   The  y?  represents  the  itn 
component  of  the  y  vector  in  the  optimal  solution  as  pro- 
vided for  in  the  Kuhn-Tucker  Theorem.   Letting  X    =    2X '  -  2X1' 
and  u  =  y'  -  y'',  the  following  results  can  readily  be  de- 
duced from  the  previous  equation  and  the  original   con- 
straints : 


X    y?  +  y 


n 


X    y°  +  y 


v 


A  y?  +  y  +  it. 


if  a.  <  x.  <  b. 

iii 


or  if  x.  =  a.  and  y:  >  0 
i    l      l 


or  if  x-  =  b.  and  y:  <  0 
i     l      '  l 


if  x.  =  a.  and  y?  =  0 
i     l      '  l 


if  x.  =  b.  and  y?  =  0 


Since  v.  >  0  and  tt  .  >  0  from  the  Kuhn-Tucker  Theorem, 
l  -        l  - 

it  is  equivalently  true  that 
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Q.    <    X 
l    - 

*i 

+  y 

Q.    >    X 
1    - 

*i 

+  y 

ft.    =    X 

1 

y\ 

+  y 

if  x.  =  b.  and  y?  =  0     (1) 
li      l        v  J 


if  x.  =  a.  and  y?  *  0     (2) 
11      l        v  ' 

otherwise  conditions  satis- 
fying the  original  con-   (3) 
straints . 


Therefore,  if  the  optimal  solution  exists,  then  X  and 
y  exist  and  are  non-negative  and  satisfy  the  above  condi- 
tions . 

The  next  step  is  to  consider  the  sufficiency  of  the 
conditions.   It  is  asserted  that  if  X  and  y  exist  and  X  is 
non-negative  and  if  y°  satisfies  results  in  (1)  ,  (2)  ,.  and 
(3)  above,  then  y°  maximizes  Eft*y. 

To  prove  this  let  y  be  any  arbitrary  direction  vector 
satisfying  the  conditions  derived  via  the  Kuhn-Tucker 
Theorem.   Observe  that  ft«y  =  Efi.Y-  =   (ft-  ~  y)y-  since 
Ey.  =  0  from  the  original  problem  constraints.   This  is  true 
for  all  i  representing  elements  of  the  ft  and  y  vectors. 
Therefore,  the  summation  can  be  broken  into  three  groups 
of  terms  as  indicated  in  (1),  (2),  and  (3)  above;  i.e., 

ft.y  =  EI(fti-y)yi  +  ZII(fti-y)y.  +  EIII(fti-y)yi 

where  Z  ,S   ,  and  E     are  summations  over  terms  in  (1),  (2), 
and  (3)  respectively.   E   <  0  since  (ft.-y)  <  0  from  the  K-T 
necessary  condition  (1)  and  y.  >  0  from  the  original  con- 
straint.  By  similar  reasoning  it  can  be  seen  that  E    <  0. 
Therefore,  it  is  true  that  ft*y  <  E 

However,  EIII(ft.-y)y.  =  XEIXIy?y.  from  K-T  condition 
'      v  l  J   !  i         '  i  '  i 

(3)  and  XEIITy?y.  =  X  I.    y?y.  since  y?  =  0  for  all  i  not 
1 i ' i     all  '  i  '  i       '  i 

1    12 


contained  in  Z    ,  and  A   ,,  Y-Y-  <  ^(  11  Y02  )  C  t-i  V2  -»    a 

*       all  'i'i  -   vall   -   ' vall  '.  )  =  X 

i  i    x     i   x 

from  the  Schwartz  inequality  and  from  the  original  con- 
straints.  Therefore,  n--y  <  A,  but  ft*Y°  =  £(G..-U-)Y?  =  AZy?2 
=  X.   And  since  A  >  0,  it  is  true  that  Q»y  <  fi~Y°  with  the 
resultant  fact  that  y°  maximizes  ff2* y  as  was  to  be  shown. 
Furthermore,  the  maximum  value  is  A  or  in  terms  of  the  suf- 
ficiency conditions;  if  A  can  be  found  then  the  optimal 
solution  is  obtained. 

At  this  point  it  is  convenient  to  rewrite  equation 
(3)  and  the  conditions  necessary  for  the  equality:: 

n.  =  Ay-  +  u  if  a.  <  x.  <  b. 

1     1  111 

or  if  x.  =  a.  and  y-  >  0 

1     1      ■  1 

or  if  x-  =  b-  and  y°.    <  0. 
1     1      '  1 

Letting  E  '  indicate  a  summation  over  the  terms  ful- 
filling these  conditions  and  summing,  results  in 

E'ft.  =  AZ'y?  +  Z'y  =  E'y 
1       '  1 

or  (4) 

where  N'  is  the  total  number  of  terms  fulfilling  the  condi- 
tions.  If  the  equation  is  squared  first  and  then  summed, 
the  expression 

A  =  E'(Q.  -y)2  (5) 

results.   It  also  follows  from  substitution  that 


Y?  ■  (^  "  lO/X.  (6) 
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With  a  method  of  calculating  values  for  u,  A,  and  y? 

now  in  hand,  development  of  the  algorithm  can  proceed. 

For  notational  ease,  let  DS  =  the  "Distinguished  Set" 

and  be  defined  as  that  set  of  indices  i  such  that  a.<x.<b., 

1   1   i* 

or  x.  =  a.  and  Y-  >  0>  or  x.  =  b.  and  y?  <    0.   In  other 
l     l      l  l     l      l 

words,  the  distinguished  set  represents  those  elements  of 
the  allocation  vector,  X,  where  the  element  does  not  lie  on 
the  boundary;  or  if  the  element  is  on  the  boundary,  then 
the  corresponding  direction  vector  element  implies  a  move 
away  from  the  boundary.   With  this  notation  equation  (4) 
can  be  rewritten  as 

where  N^-  is  the  total  number  of  elements  in  DS .. 

The  elements  contained  in  DS  can  be  separated  into 
three  categories  as  follows: 

Category  #1.   The  X  element  is  not  on  its  boundary  and 
there  is  no  restriction  on  the  y  element. 

Category  #2.  The  X  element  is  on  its  lower  boundary 
and  the  y  element  is  positive.   From  equation  (6)  then 

Q.  >  u. 

l 

Category  #3.  The  X  element  is  on  its  upper  boundary 
and  the  y  element  is  negative.   From  equation  (6)  then 

Q.    <    \L. 

l 

It  is  evident  now  that  the  algorithm  for  finding  y° , 
or  the  direction  of  fastest  increase,  must  search  out  and 
find  each  of  the  elements  included  in  the  three  categories 
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above.   Recall  that  y  is  actually  the  average  of  all  ft.  for 
which  Y-  f-    0-   This  means  it  is  necessary  to  determine  y 
concurrently  with  determining  the  elements  in  DS .   The  fol- 
lowing procedure  provides  the  method  for  doing  just  this. 

A.    THE  DIRECTION  FINDING  ALGORITHM 

1     Z 
Step  #1.   Let  y  =  ^ —   .  nc;  ft.  and  DS  =  {(f)}  to  start. 

Step  #2.   Place  into  DS  each  index  i  for  which  a.<x.<h.. 
r  ill 

Note  that  this  implies  no  restriction  on  y..   Note  also  that 

throughout  the  algorithm,  once  i  has  been  entered  into  DS , 

it  is  no  longer  considered  as  a  candidate  for  entry. 

Step  #3.   If  DS  =  {(})},  then  temporarily  let  y  =  min{ft.}. 

i 

Do  not  place  this  i  into  DS  (the  minimum  gradient  element 
serves  simply  as  a  starting  point).   If  DS  f    {<{>},  then  cal- 
culate y  from  equation  (7). 

Step  #4.   Search  indices  for  cases  where  x.  =  b..   Find 
r  11 

the  smallest  ft.  associated  with  this  search.   If  this  ft .    <  y, 
l  l 

then  enter  this  i  into  DS  and  recalculate  y.   If  this  ft   >  y, 
then  proceed  to  Step  #6. 

Step  #5.   Repeat  Step  #4  until  no  additional  i's  can 
be  entered  into  DS . 

Step  #6.   Search  indices  for  cases  where  x-  =  a. .   Find 

the  largest  ft.  associated  with  this  search.   If  this  ft.  >  y, 
&     i  l 

then  enter  this  i  into  DS  and  recalculate  y.   If  this  ft.  <  y, 
then  proceed  to  Step  #8. 

Step  #7.   Repeat  Step  #6  until  no  additional  i's  can 
be  entered  into  DS . 
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Step  #8.   Repeat  Step  #4  through  Step  #7  until  no  ad- 
ditional i ' s  can  be  entered  into  DS . 

Step  #9.   Calculate  X  -  [  I    (fi.  -  u)2J^.   If  X  =   0 

ieDS   x 
then  STOP. 

Step  #10.   Calculate  y°  from  equation  (6). 

This  concludes  the  algorithm.   At  this  point  two 
special  case  considerations  should  be  discussed: 

CASE  I.   The  possibility  exists  that  only  one  element 
will  be  entered  into  the  DS .   In  this  event  X  =  0  and  the 
optimal  solution  lies  on  the  boundary.   y°  at  this  point 
indicates  that  any  better  solution  would  violate  the  orig- 
inal constraints. 

CASE  II.   If  more  than  one  element  has  been  entered 

into  DS  and  X  =  0.  then  Q.  =  u  for  all  !►   In  this  event 

'        l 

all  ft.  are  obviously  equal  and  the  optimal  solution  has 
been  reached.   In  fact,  the  value  of  X  is  a  measure  of  how 
close  the  procedure  has  come  to  the  optimal  solution.   It 
may  be  true  in  a  practical  application  of  the  algorithm  that 
X  will  diminish  more  slowly  with  each  successive  computation 
as  the  optimal  solution  in  neared.   It  may  be  necessary  in 
this  case  for  the  user  to  make  a  determination  of  the  pre- 
cision desired  and  to  terminate  calculations  when  a  satis- 
factory degree  of  precision  has  been  reached. 

In  the  next  section  an  example  application  of  the  algo- 
rithm in  a  practically  oriented  problem  is  provided. 
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III.   EXAMPLE  APPLICATION  OF  DIRECTION  FINDING  ALGORITHM 

Consider  a  military  oriented  problem  in  which  an  ord- 
nance delivering  organization  has  responsibility  for  and 
the  capability  of  firing  its  ordnance  into  a  particular 
sector  of  territory.   Assume  that  this  organization  has 
some  means  of  assigning  probabilities  to  the  potential 
existence  of  a  target  at  any  of  a  number  of  specific  loca- 
tions within  its  sector  of  responsibility.   Assume  also 
that  the  organization  knows  what  damage  it  can  expect  from 
a  specific  type  round  delivered  at  a  specific  point.   Know- 
ing this  information,  the  organization  would  then  like  to 
know  how  to  optimally  distribute  the  inventory  it  intends 
to  deliver  so  that  the  probability  of  target  survival  is 
minimized. l 

As  a  first  step  to  solution  of  this  problem,  the  sector 
of  responsibility  is  assumed  to  be  square  in  geometric  shape 
and  an  imaginary  grid  is  imposed  upon  the  sector.   The  in- 
dividual squares  of  the  grid  are  numbered  sequentially  be- 
ginning in  the  top  left  corner  and  working  first  across  and 
then  down.   Location  of  a  target  (or  the  possible  location 
of  a  target)  can  then  be  referenced  by  a  specific  numbered 
square  of  the  grid.   Expected  damage  in  any  grid  square  can 
be  calculated  from  the  number  of  rounds  which  impact  in  that 
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This  problem  is  formulated  and  presented  for  solution 
in  a  thesis  paper  prepared  by  Captain  William  A.  Hesser, 
USMC,  at  the  U.  S.  Naval  Postgraduate  School  in  February/ 
March  1971. 
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square  and/or  in  nearby  squares.   Allocation  of  ordnance 
is  represented  by  the  number  of  rounds  delivered  to  a 
specific  grid  square. 

The  survival  probability  of  a  target,  or  the  objective 
function  which  is  to  be  minimized,  can  now  be  expressed  as 
a  mathematical  equation 


f(y)  =  z  p.exp{-z  e^y,} 

i        J    J  J 


where  the  variables  have  interpretations  as  follows: 

p.   =  probability  that  target  is  located  in  it A  square, 

$..  =  probability  of  destruction  of  a  target  appearing 
in  the  i    square  by  a  round  detonating  in  the 
jtn  square  (this  is  referred  to  as  the  damage 
function)  , 
y.   =  percent  of  total  inventory  to  be  expended  allo- 
cated to  the  j tn  square, 
F(Y)  =  survival  probability  of  the  target  and  is  a 
function  of  the  allocation  of  Y. 
Note  at  this  point  that  the  survival  function,  F(Y), 
is  a  convex  function  and  the  variable  Y  is  actually  a 
closed  convex  set  (there  are  only  a  finite  number  of  rounds 
which  can  be  fired) .   Also  note  that  in  the  trivial  case 
when  Y  =  (0,0,...,0),  i.e.,  no  firing  takes  place,  the 
survival  function  reduces  to  F (Y)  =  .  p.  =  I  which  is  ap- 
propriate . 

The  problem  now  is  to  determine  what  allocation  of  Y 
will  minimize  F(Y)  . 
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It  can  be  seen  at  this  point  that  in  order  to  obtain 
realistic  and  usable  results,  it  is  necessary  to  use  a  grid 
with  a  large  number  of  squares.   Correspondingly,,  the  allo- 
cation vector  is  composed  of  a  large  number  of  elements. 

The  problem  can  now  be  represented  as  follows: 
minimize  F (Y)  =  E  p . exp {-£$  .  .y . } 

subject  to  constraints 

I  y.  =  1.0 

.  '  i 

l 

and 

y.  >  0.0. 
3  i  - 

From  previous  discussion  and  assuming  the  p.  and  $- 

are  known  beforehand,  this  problem  can  be  solved  using  the 

Direction  Finding  Algorithm.   Observe  that  the  first  partial 

derivative  required  in  the  use  of  the  algorithm  is 


F   (Y)  =  Z  p.  6- •  exp{-Z  g. .y. }. 


Solution  of  this  problem  via  the  Direction  Finding  Algo- 
rithm was  accomplished  on  the  IBM  360  computer  system  at  the 
U.  S.  Naval  Postgraduate  School.   Programming  was  done  in 
FORTRAN  IV  Language  using  the  G-Level  Compiler.   A  copy 
of  the  program  listing  is  included  after  the  appendices. 

A.    EXPLANATION  OF  EXAMPLE 

Since  the  purpose  of  this  example  is  to  demonstrate  ap- 
plication of  the  algorithm  and  not  to  solve  a  given  existing 
problem,  assignment  of  values  to  parameters  of  the  problem 
was  based  on  academic  interest  and  convenience  rather  than 

practically  oriented. 
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A  grid  size  of  3600  squares,  or  equivalently  a  60X60 
matrix  was  used.   It  is  felt  this  grid  size  selection  is 
consistent  with  what  might  be  used  in  a  truly  practical 
application  and  at  the  same  time  it  enlarges  the  problem 
sufficiently  so  that  information  about  machine  time  for 
calculation  would  be  meaningful. 

The  distribution  of  probabilities  (p.'s)  used  in  this 
example  was  uniform  in  nature.   A  p  =  1/3600  was  assigned 
to  each  of  the  3600  grid  squares.   In  a  normal  application 
it  is  expected  that  most  grid  squares  would  have  a  p  =  0.0 
assigned  to  them  and  only  a  relatively  small  number  would 
be  assigned  a  positive  probability.   Those  positive  assign- 
ments would  be  most  likely  based  on  observer  reports,  known 
routes  of  movements,  specific  areas  offering  good  cover  and 
concealment,  and  other  considerations  such  as  these.   How- 
ever, in  this  example  a  uniform  distribution  was  selected 
for  two  reasons.   Firstly,  making  all  p's  positive  is  con- 
sidered a  worst  case  condition  with  respect  to  the  number 
of  individual  calculations  required  and  the  resultant  ma- 
chine time  necessary  for  execution.   Since  solution  by  the 
Direction  Finding  Algorithm  is  an  iterative  process,  machine 
time  is  an  important  consideration  and  a  known  upper  limit 
on  this  factor  is  extremely  meaningful.   Secondly,  a  uniform 
distribution  was  used  so  that  effects  on  the  boundaries  of 
the  grid  could  be  observed.   Since  the  damage  function  (3--'s) 
implies  interaction  between  grid  squares,  it  is  not  immedi- 
ately obvious  what  the  optimal  allocation  along  the  edges  of 

the  grid  should  be. 
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Selection  of  an  appropriate  damage  function  would  be 
based  on  the  type  of  weapon  being  employed  and  the  result 
of  statistical  data  obtained  for  that  weapon..   Inherently, 
the  damage  function  selected  directly  affects  total  machine 
time  necessary  to  solve  the  problem.   However,  this  is  a 
variable  which  is  dependent  on  a  specific  situation  and 
not  of  prime  consideration  in  this  paper*   A  simple  expon- 
ential relationship  was  selected  for  the  example  because 
it  was  felt  that  an  exponential  function  would  be  repre- 
sentative of  most  damage  functions  in  terms  of  accuracy 
and  machine  time  necessary  for  calculation. 

Given  the  above  explanation  of  problem  parameters  it 
now  becomes  necessary  to  develop  a   computer  oriented  logi- 
cal approach  to  solution  of  the  problem. 

The  logical  sequence  followed  in  finding  the  direction 
of  fastest  increase  is  explained  earlier  in  this  paper.   At 
this  point  it  should  be  noted  that  the  program  actually  em- 
ploys two  Direction  Finding  Algorithms.   In  one  case 
(SUBROUTINE  DF$ALL   in  the  program  listing)  the  y   or  direc- 
tion element  associated  with  each  and  every  element  of  the 
allocation  vector  is  calculated.   This  is  necessary  before 
a  final  optimal  solution  can  be  obtained.   However,  while 
working  toward  the  optimal  solution,  a  significant  reduction 
in  calculation  time  can  be  realized  by  finding  the  direction 
of  fastest  increase  only  on  the  variables  which  have  a  posi- 
tive allocation  currently  associated  with  them  or  a  positive 
Y  in  case  the  allocation  is  currently  zero.   For  this  reason 


21 


the  program  maintains  a  list  (SUBROUTINE  LIST)  of  those  ele- 
ment indices  which  meet  this  criteria.   With  each  iteration 
then,  SUBROUTINE  DF$LST  is  employed  to  find  the  direction 
of  fastest  increase  only  for  those  elements  appearing  on 
the  "list."   When  the  optimal  solution  is  found  for  the 
current  "list,"  execution  of  the  program  shifts  to  SUB- 
ROUTINE DF$ALL.   After  execution  of  DF$ALL ,  a  new  list  is 
formed  and  the  above  outlined  process  is  repeated..   Only 
when  optimization  is  indicated  through  the  use  of  DF$ALL 
is  execution  terminated. 

In  order  to  employ  the  Direction  Finding  Algorithm  as 
indicated,  various  other  calculations  obviously  must  be 
made.   Other  routines  in  the  program  are  written  to  accom- 
plish tasks  such  as  evaluating  the  gradient,  evaluating 
the  objective  function,  and  revising  the  current  allocation 
to  one  which  is  closer  to  the  optimal  solution*   These  tasks 
must  necessarily  be  carried  out  during  each  iteration  of 
the  solution  process.   A  general  flow  chart  outlining  the 
logical  procedure  followed  is  depicted  in  Appendix  A.   A 
complete  listing  of  subroutines  employed  along  with  a  des- 
cription of  the  purpose  of  each  is  provided  in  Appendix  B. 

Results  obtained  in  the  solution  to  this  example  prob- 
lem are  included.   It  is  interesting  to  note  in  these  results 
that  the  optimal  solution  has  no  allocation  of  resources  in 
the  rows  or  columns  immediately  adjacent  to  the  edges  of  the 
grid.   Obviously  the  damage  resulting  in  these  grid  squares 
from  rounds  allocated  to  and  impacting  in  nearby  interior 


22 


grid  squares  is  considered  sufficient  to  compensate  for  loss  of 
damage  effect  that  would  occur  if  rounds  were  allocated  to 
the  outermost  rows  and/or  columns.   The  higher  allocation  in 
those  nearby  interior  grid  squares  suggests  this  conclusion. 

It  is  also  interesting  to  note  that  the  optimal  solu- 
tion is  not  completely  symmetric  as  might  be  expected  from 
using  a  uniform  type  distribution  of  the  p.  *  s ..   Inspection 
of  the  final  allocation  shows  the  upper  left  corner  of  the 
array  somewhat  different  from  allocations  assigned  in  the 
other  three  corners.   This  is  the  result  of  the  starting  al- 
location used  in  the  solution  process  and  the  precision  con- 
sidered acceptable  within  the  program.   The  initial  starting 
feasible  solution  was  for  the  entire  inventory  to  be  assigned 
to  the  extreme  upper  left  grid  square.   This  heavy  allocation 
in  just  one  location  prevented  a  lesser  positive  allocation 
from  appearing  in  nearby  grid  squares  as  the  program  stepped 
along  towards  optimality.   The  program  was  actually  moving 
toward  a  totally  symmetric  final  allocation,  but  settled  for 
something  less  because  of  the  precision  cutoff  programmed 
into  the  problem.   If  a  higher  degree  of  precision  were  de- 
manded and  if  the  program  were  rewritten  to  provide  greater 
machine  accuracy,  the  program  would  tend  to  a  completely 
symmetric  optimal  solution. 

Certain  difficulties  were  encountered  in  programming 
this  problem  for  machine  solution.   A  discussion  of  some  of 
these  difficulties  and  comments  on  observations  is  considered 
of  academic  interest  to  the  reader  at  this  point. 
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Since  use  of  the  Direction  Finding  Algorithm  is  an 
iterative  process  which  moves  from  any  feasible  starting 
solution  to  a  better  one  and  then  ultimately  to  an  approxi- 
mate optimal  solution,  selection  of  an  initial  starting 
point  can  be  a  critical  matter.   Initial  attempts  at  solving 
this  problem  were  done  by  assigning  the  entire  inventory 
to  the  first  grid  square  and  proceeding  from  there.   Anal- 
ysis of  results  obtained  from  this  "start  from  scratch"  ap- 
proach indicated  that  machine  time  in  excess  of  twelve  hours 
would  be  required  to  completely  solve  the  problem.   This  re- 
quired time  was  considered  unacceptable  and,  therefore,  an 
alternate  approach  was  devised.   The  sector  covered  with  a 
3600  grid  was  initially  considered  to  be  covered  with  a  400 
element  grid.   An  optimal  solution  was  found  for  the  400  grid 
case  and  this  solution  was  then  used  as  a  starting  solution 
for  the  3600  grid  case.   Results  in  terms  of  machine  time 
were  remarkable.   The  optimal  solution  to  the  400  case  was 
obtained  in  17.5  seconds  of  machine   execution  time.   Over- 
all machine  time  to  the  final  optimal  solution  was  six 
minutes  nine  seconds . 

A  point  of  interest  in  the  program  is  the  value  of  "D" 
or  the  distance  which  the  allocation  is  moved  in  each  itera- 
tion.  Recall  the  lower  and  upper  boundaries  of  the  alloca- 
tion elements  are  0.0  and  1.0  respectively.   The  lower  bound 
on  D  is  0.0  which  is  synonomous  with  no  move  at  all.   The 
upper  bound  on  D  is  the  maximum  distance  across  the  Y-space 
and  is  /~~2 .   In  some  calculations  within  the  program  D  is 
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controlled  by  a  variable  which  is  near  its  boundary  and 
trying  to  move  in  the  direction  of  its  boundary..   The 
result  of  this  situation  in  large  problems  such  as  this 
example  is  that  D  becomes  very  small  and  this  leads  to 
short  moves  on  each  iteration.   When  many  variables  are 
close  to  their  boundaries  and  trying  to  move  in  that  direc- 
tion, the  overall  impact  is  a  significant  increase  in 
machine  execution  time.   However,  there  are  various  places 
in  the  solution  process  where  D  can  be  successfully  reset 
to  a  predetermined  value  so  that  time  is  not  wasted  on 
relatively  small  moves  when  larger  moves  are  permitted. 
There  is  also  the  problem  of  what  to  do  with  D  when  move- 
ment is  made  in  a  direction  of  increase  but  it  has  not  re- 
sulted in  an  improved  value  of  the  objective  function.   To 
resolve  these  difficulties,  an  initial  D  value  of  0.5  was 
selected.   This  had  the  effect  of  making  the  first  step  a 
large  one.   A  reset  value  of  0.1  was  decided  upon  and  each 
time  a  new  list  was  compiled,  D  was  reset  to  this  value. 
A  modification  procedure  was  used  to  adjust  D  when  a  move 
was  made  but  an  improvement  in  the  objective  function  was 
not  realized  (rightly  enough  referred  to  as  a  "failure") . 
In  this  case  D  was  simply  cut  in  half  and  the  move  tried 
again.   In  the  case  of  success  in  improving  the  value  of  the 
objective  function,  D  was  left  alone  and  the  same  value  used 
in  the  next  iteration.   The  only  exception  to  this  procedure 
was  in  the  case  of  a  "success"  which  had  been  immediately 
preceded  by  a  failure.   In  this  circumstance  D  was  halved 
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for  the  next  iteration.   To  understand  the  rationale  be- 
hind this  consider  the  following  simplified  two-dimensional 
representation  of  a  concave  function  which  is  to  be  maxi- 
mized: 


Let  point  A  represent  the  current  allocation  and  the  direc- 
tion of  fastest  increase  is  obviously  to  the  right.   Sup- 
pose the  present  D  is  equal  in  value  to  (B-A)  ..  When  the 
move  to  B  is  attempted,  a  "failure"  occurs  since  the  value 
of  the  function  is  less  than  it  was  at  point  A.   The  pro- 
cedure now  calls  for  halving  D  and  try  to  move  again.   This 
time  the  move  is  to  point  C  and  a  "success"  is  realized. 
If  D  is  not  adjusted  at  this  point,  the  next  iteration 
would  call  for  an  attempted  move  back  to  point  A  which  was 
the  starting  point  for  the  previous  iteration.   To  preclude 
this  return  to  a  previous  allocation,  D  is  cut  in  half 
after  the  occurrence  of  a  success  which  had  been  immedi- 
ately preceded  by  a  failure.   It  should  be  noted  that  this 
procedure  becomes  particularly  useful  when  the  optimal 
solution  is  being  neared  and  oscillations  around  the  optimal 
solution  start  occuring. 
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Another  consideration  worthy  of  note  is  the  influence 
of  machine  accuracy  on  the  results  of  the  problem.   In  a 
problem  of  this  sort  with  potentially  3600  variables  posi- 
tive simultaneously,  the  effects  of  calculating,  founding, 
and  cumulatively  adding  these  variables  must  be  taken  into 
account.   Early  attempts  at  solution  were  made  using  single 
precision  under  FORTRAN  IV.   This  approach  was  quickly 
changed  to  double  precision  for  all  real  variables  used  in 

the  program.   Even  so,  results  obtained  were  considered 

_  i  o 
suspect  in  accuracy  beyond  10    .   To  illustrate  the  prob- 
lem, in  this  example  the  initial  p.  value  assignment  to 
each  grid  square  was  made  by  assigning  each  grid  square  a 
value  of  1/400.   This  value  is  supposed  to  be  accurate  in 
FORTRAN  IV  to  10      As  a  check  device  the  program  then 
sums  all  p's  assigned  and  prints  the  results.   The  re- 
sultant sum  of  400  separate  additions  differs  from  1.0  by 
2.24x10   .   Whether  or  not  accuracy  at  this  level  is  tol- 
erable is  up  to  the  user,  but  most  certainly  is  worthy  of 
some  consideration. 

Machine  storage  space  and  execution  time  requirements 
in  this  problem  were  sufficiently  significant  to  be  worthy 
of  comment.   Because  of  the  requirement  to  maintain  a  large 
number  of  large  arrays  of  stored  information  (in  the  common 
storage  region  alone  there  are  six  arrays  each  of  which  con- 
tains 3600  double  precision  real  numbers) ,  this  example  pro- 
gram run  on  the  IBM  360  System  required  352',  000  bytes  of 
storage.   Some  space  saving  techniques  were  employed  and, 


27 


undoubtedly,  additional  ones  could  be  found.   However,  it 
should  be  recognized  that  the  storage  requirements  are  not 
a  direct  result  of  the  methodology  used  in  the  Direction 
Finding  Algorithm,  but  rather  the  result  of  the  nature  of 
the  problem  being  solved.   As  mentioned  previously,  total 
execution  time  necessary  was  six  minutes  nine  seconds.   The 
approach  of  solving  the  problem  on  a  small  scale  first  and 
then  using  these  results  as  a  starting  solution  for  the 
larger  case  was  an  instrumental  factor  in  keeping  execution 
time  within  reason.   Other  techniques  involving  the  applica- 
tion of  "clean  coding"  contributed  toward  reducing  time 
requirements.   One  technique  developed  is  considered  to  be 
of  academic  interest  and  involved  the  calculation  of  geo- 
metrical distance  between  squares.   This  calculation  is 
necessary  in  evaluation  of  the  gradient  and  in  evaluating 
the  function  value.   The  distance  calculation  involves  sim- 
ply measuring  the  horizontal  and  vertical  distance  between 
the  two  grid  squares  in  question  and  then  evaluating  the 
square  root  of  the  sum  of  squares.   To  illustrate  the  signi- 
ficance of  this  simple  routine,  it  is  done  in  excess  of  one 
million  times  in  each  iteration  of  the  program.   The  execu- 
tion time  required  to  call  the  built-in  machine  square  root 
function  and  evaluate  the  argument  is  known  to  be  in  excess 
of  200  microseconds.   Since  this  obviously  contributes  to 
exorbitant  machine  execution  time,  it  was  decided  to  make 
the  distance  computation  between  any  two  grid  squares  once 
and  for  all  at  the  beginning  of  the  program  and  then  store 
this  information  for  later  reference  as  required.   This  was 
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possible  since  the  distance  evaluations  were  somewhat  re- 
petitive in  nature.   The  exact  time  saving  incurred  by  ap- 
plication of  this  technique  was  not  evaluated.   However,  a 
conservative  estimate  is  that  execution  time  was  reduced 
by  a  factor  of  five  and,  therefore,  makes  the  nominal  ad- 
ditional storage  requirement  well  worth  while.   The  specific 
operation  just  described  is  contained  in  SUBROUTINE  REFER. 
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IV.   SUMMARY  AND  CONCLUSIONS 

The  Direction  Finding  Algorithm  can  be  applied  to  a 
variety  of  concave  or  convex  dif ferentiable  functions  to 
find  an  optimal  solution.   However,  the  domain  of  these 
functions  must  be  represented  as  a  closed  convex,  set. 
Its  applicability  in  practical  problems  is  particularly 
useful  when  the  objective  is  to  optimize  allocation  of 
resources . 

In  large  scale  problems  it  is  easily  adaptable  to 
machine  utilization.   Time  and  storage  requirements  for  a 
machine  solution  are  primarily  dependent  on  the  nature  of 
the  specific  problem  being  solved.   The  algorithm  has 
relatively  little  effect  on  machine  storage  requirements. 
However,  the  algorithm's  impact  on  machine  time  is  the  re- 
sult of  the  boundary  problem  previously  discussed.   A 
method  of  overcoming  the  boundary  problem  as  well  as  an 
extended  application  of  this  algorithm  into  an  area  of 
problems  which  require  the  allocated  variable  to  be  repre- 
sented in  typical  matrix  form  rather  than  in  vector  form 
is  currently  being  investigated  by  the  aforementioned 
Dr.  Danskin. 
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APPENDIX  A:   FLOWCHART -GENERAL  LOGICAL  PROCEDURE  FOLLOWED 


ADJUST 
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APPENDIX  B:   SUBROUTINE  LISTING 

SUBROUTINE  AMMCAL :   AMMCAL  evaluates  the  objective  function. 
SUBROUTINE  AREA:   AREA  determines  which  jrs  are  associated 
with  a  given  i  in  the  3- •  term. 

SUBROUTINE  CHANGE:   CHANGE  converts  the  optimal  solution 
obtained  from  the  400  grid  case  into  the  starting  solution 
for  the  3600  grid  case.   It  also  changes  values  of  various 
parameters  used  within  the  program  for  use  in  the  3600  grid 
solution. 

SUBROUTINE  CNTROL :   CNTROL  controls  the  sequencing  of  events 
in  solution  of  the  problem.   It  tests  each  attempted  move 
for  success  or  failure  and  adjusts  distance  (D)  accordingly. 
SUBROUTINE  DF$ALL:   DF$ALL  is  the  Direction  Finding  Algo- 
rithm which  finds  the  direction  of  fastest  increase  for  the 
entire  allocation  vector. 

SUBROUTINE  DF$LST:   DF$LST  is  the  Direction  Finding  Algo- 
rithm which  finds  the  direction  of  fastest  increase  just 
for  those  variables  on  the  list. 

SUBROUTINE  LIST:   LIST  compiles  a  list  of  those  variables 
which  are  positive  and/or  those  variables  which  have  a 
positive  y. 

SUBROUTINE  MOVE:   MOVE  adjusts  the  allocation  vector  to 
reflect  an  attempted  move  in  the  direction  of  fastest  in- 
crease.  It  moves  to  the  boundary  any  variables  which  are 
within  10     of  their  boundary. 
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SUBROUTINE  OMCALC :   OMCALC  calculates  the  vector  of  first 
partial  derivatives  of  the  allocation  vector. 
SUBROUTINE  REFER:   REFER  fills  a  matrix  with  values  rep- 
resenting geometrical  distance  between  two  grid  squares. 
SUBROUTINE  SET$P:   SET$P  sets  precision  desired  into  the 
program  and  fills  the  P  vector. 
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HEPE  TS  THF  MAIN  PROGRAM 

YY  IS  THF  3ASE  POINT 

J  IS  A  POINTER  TO  ELEMENTS  ON  THE  LIST 

LONG  IS  THF  TOTAL  NJMRER  OF  ELEMENTS  IN  VECTORS 

SI7F  IS  LENGTH  OF  OME  SIDE  OF  SQUARE  G«?Tn 

FAIL1  IS  TRUE  IF  SUCCESS  FOLLOWS  A  FAILURE 

SHORTD  IS  TRUE  IF  WORKING  WITH  A  D  CLOSE  TO  THF  BOUNDARY 

JJJ  IS  1  IF  ABBREVIATED  OMCALC  IS  DESIRED 

IMPLICIT  RFAL*9( A-H,n-z,$) 

LOGICAL  SHORTDtFAILl 

COMMON  Y(3  600) , YY( 3 600) ? J ( 3600) ,GAf 36001 rnM( 3600) , 
1 IMC 3600 ),P( 3630 )  ,E<  36  00)  T  DSTNCEC  8  ,  8  )  , 0, V? , AMM , AM, KO, 
1L0NS,III  ,JJJtJK,TSIZ<=,  I  RANGE,  ILO,  IMI  T  IT  ED  ,  SHORTO,  FA  IL  1 
BLANK  COMMON  STATEMENTS  IN  SUBROUTINES  HAVE  8EEN  OMITTED 
FOR  THIS  PRINTING 

11=0 

111  =  0 

JJJ  =  0 

KKK  =  0 

KKKK=0 

NN  =  0 

ITER=0 

LONG=400 

ISI?F=20 

IPANGF=2 

IL0=3 

IHI=18 

AM=1.0 

0=0.1 

FAIL1=. FALSE. 

SHORTD=. FALSE. 

DO  1  1=1,3600 

Y(I)=0.0 

YY(I )=0.0 

OM(T )=n.o 

GA( I )=0.0 

TM( T)=n 

J(I)=0 
1     CONTIMUF 
FOLLOWING  TWO  CARDS  INDICATE  STARTING  SOLUTION 

Y(l)=1.0 

YY(1)=1.0 

CALL  SCT4;P 

CALL  REFER 

CALL  CNTROL 

CAIL  CHANGE 

CALL  CNTPOL 
AT  THIS  oniNT  PROBLEM  TS  SOLVFD  E  NECESSARY  WRITE 
STATEMENTS  NFED  TO  BE  INSERTED 

STOP 

END 

SUBRHUTIMF  ONTROL 
THIS  PORTION  CONTROLS  EXECUTION  OF  THE  PROGRAM 

IMPLICIT  RPAL*8( A-H,0-Z,$) 

LOGICAL  SH0PTD,FAIL1 
1000  CALL  LIST 

CALL  AVMf.AL 
FIRST  timp  THROUGH  AMMCAL  TT  IS  CUT  off  SHORT 

CALL  OMCALC 

CALL  OE<;ALL(Y,OM,GA,LONG,V2,r.lO<?0) 
AT  THIS  POINT  GET  0^  WITH  NORMAL  EXECUTION 
1010  CALL  LIST 

0=0.1 

IF(III.FO.O)  5=0.5 

111  =  1 

FAIL1=. FALSE. 

SHORTD=.FALSF. 
1015  DO  1020  K=1,K0 

JK=J(K) 

YYY=YY( JK) 
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GAA=GAUK) 

IF(YYY+D*G4A.GE.0.0)  GO  TO  1020 

D=-YYY/GAA 

SHOPTp=.TRUE. 
1020  CPNTINUP 
1025  CALL  MOVE 

CALL  AMMCAL 

TF(AMM.LT.AM)  50  TD  1040 
BRANCH  TO  104H  BELOW  IP  SUCCESS-OTHEPVfl  SE  DO  FOLLOWTNG: 

0=0/2 

FAIL1=.TRUP. 

GO  TO  1025 
IF  SUCCESS, PT°ST  SET  NEW  AM  E  MOVE  BASE  POINT 
1040  AM=AMM 

DO  1045  K=1,K0 

JK=J(K) 
1045  YY(JK)=Y(JK) 

IF(FAIL1 )  GO  TO  1050 

GO  TO  1055 
1050  n=0/2 

FAIL1=. FALSE. 
10^5  GALL  OMCALC 

CALL  DF$LST(v,0MtpA,J,KQ,LONGtV2,E1060) 

IF(.NOT.SHORTD)  GO  TD  1015 

0=0.1 

SHORTD=. FALSE. 

GO  TO  1015 
1060  JJJ=1 
HFPE  CALCULATE  ONLY  THOSE  0MEGAS  NOT  ALREADY  CALCULATED 

DO  1070  JJ=1,L0NG 

IF(  IM( JJ).EQ.l  )  GO  TO  1070 

OM(JJ)=0.0 

JK  =  JJ 

CALL  OMCALC 
1070  CONTINUE 

JJJ  =  0 

CALL  DF $ ALL (Y,PM,GA,L0NG,V2, £1090) 

GO  TO  1010 
1090  PPTURN 

END 

SUBROUTINE  SFT$P 
THIS  PPUTTNF  EMTFPS  D  VALIFS  AT  START  OF  PROGRAM 
IMPLICIT  RCAL*8(  A-H^-j,  $) 

LOGICAL  SHHRTP^FAILl 

PREC=10.0**(-3 ) 

V=PREC/2**.5 

V2=V*V 
20    VAL=1. 0/400.0 

DO  22  N=l,400 
22    P(N)=VAL 

SUMP=0.0 

DO  25  N=l,400 
25    SUMP=SUMP+P(N) 

DO  27  N  =  1,a.00 
27    P(N)=P( N)/SUMP 

RFTURN 

END 

SUBROUTINE  LIST 
THIS  IS  LIST  CONSTRUCTING  SUBROUTINE 
IMPLICIT  °FAL*8( A-H,0-Z,$) 
LPGICAL  SHORTDtFAILI 
K=0 

ITER=ITEP+1 
DO  60  JJ=1,L0MG 
IF(YUJI.GT.O.O)  GO  TO  *0 
IF(GA( JJ) .GT.D.O)  GO  T0  50 
GO  TO  ftO 
50    K=K+1 
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J(K)=JJ 
60    CONTINUF 

KG=K 
KO  NOW  CONTAINS 
0=0.1 
RETURN 
END 


LIST  SIZE 


THIS  I 
DFSLST 
S  IS 
0M$  I 

GAS  I 
LONGS 
V2$  I 
JS  IS 
KO*  I 
*  IS 


THIS 
100 


110 
115 

120 

140 


150 


160 


SUBROUT 

S  UIREC 

ACCEPT 

VECTOR 
S  VECTO 
S  GAMMA 

IS  CUP 
S  V*V  W 

VECTOD 
S  A  MUM 
SPFCI AL 
IMPLICT 
REAL*8 
LOGICAL 
OIMENSI 
SUBPOUT 
COMMON/ 
N=0 
KK  =  0 
S  =  0.0 
DO  110 
JK=J$( K 
GA$( JK) 
I  NO  UK) 
IF(S< JK 
I  F  (  S  (  J  ;< 
INOUK) 
N  =  N  +  1 
S=s+OM<; 
CONTIN'J 
IFCN.EQ 
MU=S/N 
AOM^MU 
GO  Tn  i 
AOM=PM$ 
DO  1?0 
JK=JS(K 
OMM=OMS 
IF(OMM. 
CONTINU 
MU=AOM 
JJ  =  0 
DO  150 
JK=JS(K 
IF(IND( 
IF( $( JK 
IF(OM$( 
AOM=0MS 
JJ  =  JK 
CONTINU 
IF( AOM. 
KK=0 
MU=(N*M 
AOM=MU 
N=N+1 
IND(JJ) 
GO  TO  1 
JJ  =  0 
DO  165 
JK=JS(K 
IF(IND( 
IF($( JK 
IF(OMS( 


INE  DF$L ST rS,OM$,GA$,J$,KO$,LONG$TVZ $'»•*) 

TION  FINDING  ALGO°!THM  USING  THE  LIST 
S  AS  PARAMETERS  THE  FOLLOWING: 
OF  VARIABLE  TO  BE  ALLOCATED 
P  0^  FICST  PARTTALS 

VECTOR 
RENT  VECTOR  LENGTH 
HFOF  V=PPFCISI0N/?**.5 

OF  POINTERS  to  ELEMENTS  ON  THP  LIST 
BE°  INDICATING  THE  LIST  LENGTH 

RETURN  POINT  (RFTUPN  1} 
T  RPAL*R( A-H,0-Z,S) 
LA,MU 

I  NO 
ON  $( LONGS) T CMS (LONGS), GAS (LONGS) , J S( LONGS) 
INE  P^OUIRES  NO  BLANK  COMMON  ACCESS 
W0RKA/IND(3600) 


K=1TK0S 

) 

=  0.0 

=. FALSE. 

) .EO.n.O) 

).E0.1 .0) 

=.T^UE. 

(JK) 

F 

.0)  GO  tq 


TO 
to 


110 
110 


115 


40 

(  J$(D) 

K=1,K0$ 

) 

(  JK) 

LT.AOM) 

F 


AOM=OMM 


K=1TK0S 

) 

JK))  G3  TO 

).LT.1.0)  G 

JK) .GE. AOM) 

(JK) 


GE.MU)  GO  T 
U+A0M)/(N+1 


=.TRUF. 
40 

K=1,KC$ 

) 

JK) )  GO  TH 

).GT.O.O)  G 

JK) .LF. AOM) 


1^0 

0  TO  1^0 
GO  TO  150 


0  170 
) 


165 

0  TO  165 
GO  to  165 
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AOM=PM$( JK) 
JJ  =  JK 
165   CONTINUE 

IF ( AOM.Lc.MU)  GO  TO  ISO 

KK=0 

MU=(N*MU+AOM)/(N+l) 

AOM=MU 

N=N  +  1 

TNO(JJ)=.TRUE. 

GO  TD  160 

IE(KK.EQ.l)  GD  Tn  igo 

GO  TO  160 

KK=  1 

GO  TO  140 

SS=0.0 

00  1^5  K=ltKOS 

JK= J$ (  K  ) 

if(  ,not.tnd(  jo  )  go  to  195 

sd=pms(  jk)-mu 

ss=ss+sp*sd 

continue 

if(ss.lt.v2$)  return  1 

la=dsort(ss) 

DO  200  K=1,K0S 
JK  =  J S ( K  ) 

IF(  .NOT.INPUO  )  GO  TO  200 
SD=PMS( JK)-MU 
GAS(JK)=SD/LA 
200   CONTINUE 
RETURN 
END 


170 
180 
IPO 


195 


DF$ALL 
$  IS 
OMS  I 
GAS  I 
LONGS 
V2S  I 
*  IS 


THIS  S 

IND(I) 

500 
N  IS  C 

KK  IS 

FIRST 


510 
AT  THI 


IF  DS 
515 


SURPOUT 

ACCEPT 

VFCTPR 
S  VFCTp 
S  'GAMMA 

IS  nun 
s  v*v  w 

SDEOIAL 
IMPLICT 
RFAL*8 
LOGICAL 
DIMFNSI 
U3R0LITI 
COMMON/ 
IS  SET 
N=0 

OUNTFR 
KK  =  0 
LOGIC  S 
S  =  0.0 
TO  CHEC 
DO  510 
GA$(  J  J) 
IND( JJ) 
IF($( JJ 
IF(S( JJ 
IND< JJ) 
N=N+1 
S=S+OMS 
CONTINU 
IF(N.EO 
MU=S/N 
S  STAGP 
AOM=MU 
GO  TO  5 
IS  E^PT 
AOM=OM$ 
DO  520 
OMM=CMS 


INF  DFSALLf  *»^S*GA$,L0NS*rV?$r*1 

S  AS  PARAMETERS  THE  FOLLOWING: 
OF  TH^  VARIABLE  T0  BE  ALLOCATFD 
P  pF  FIP^T  ?A°TIALS 

vccjnt? 
RENT    VECTOR     L^NGTH 

HFRF    V=PRECISI0N/?**.5 

RCTlfPNj  POT  NT  (RETIJPN  1) 
T  REAL*8< A-H,0-Z,S) 
LA,MU 

TNP 
ON  * (LONGS) , OMS (LONG $) ,GA$( LONGS) 
NF  REQUIRES  NO  3LAMK  COMMON  ACCESS 
W0RKA/IN0C3600) 

TRUE  WHFN  $(I)  ENTERS  DS 

FOR  NUMBER  HF  ELEMENTS  IN  DS 

WITCH  WHICH  CONTROLS  FLOW 

K  FOP  FLEMENTS  NOT  ON  THEIR  BOUNDARIES 

JJ=1»L0NG$ 

=  0.0 

=. FALSE. 

) .EG. 0.0)  GO  TP  510 

) .EO.l.O)  GO  TO  510 

=.TRUE. 

(  JJ) 

F 

.0)  GO  TO  515 

,  MU  IS  AVERAGE  OF  INTERIOR  ELEMENTS 

40 

Y  SET-FIND  MINIMUM  OMFGA 

(1) 

JJ=1, LONGS 

(  JJ) 
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IMOMM.LT. AOM)  AOM  =  OMM 
520   CONTINUE 

MU=AOM 
NOW  TO  CHECK  FOR  FLFMFNTS  ON  UPPER  BOUNDARY  AND  NOT  IN  DS 
540   JMARK=0 

DO  550  JJ=1, LONGS 

IF(IND( JJ) )  GO  TO  550 

IF($( JJ) .LT.1.0)  GO  TQ  550 

If=(OM$<  JJ)  .GE.  AOM)  GO  tq  550 

AOM=OMSUJ) 

JMARK=JJ 
550   CONTINUE 

IF( AOM.GE.MU)  GO  TO  570 

KK=0 
NOW  TO  ^CALCULATE  MU 

MU=(N*^U+AOM)/f N+l) 

AOM--=MU 

N  =  N+1 

INDf JMARK)=.TRUE. 

GO  TO  540 
NOW  TO  CHECK  POP  ELEMENTS  ON  LOWER  BOUNDARY  AND  NOT  IN  DS 
560   JMARK=0 

DO  5*5  JJ=1, LONGS 

IF(IND( JJ) )  GO  TO  565 

IF(S< JJ) .GT.O.O)  GO  TQ  565 

IF(OMS( JJ) .LE.AOM)  GO  TO  565 

AOM=OMS(  JJ) 

JMAOK=JJ 
56  5   CONTINUF 

IF( APM.LE.MU)  GO  TO  SRO 

KK  =  0 
RECALCULATE    MU    AGAIN 

MU=(N*mu+A0M)/(N+1) 

AOM=MU 

N  =  N  +  1 

IND( JMARK)=.TRUE. 

GO  TO  560 
570   I^(KK.FO.l)  GO  TO  ^90 
KK=1  IMPLIES  NO  MOPE  SEARCHING  IS  PEOUrRFD 

GO  TO  560 
580   KK=1 

GO  TO  540 
NOW  TO  MND  SUM  OP  SQUARES  OP  ELEMENTS  IN  DS 
590   SS=0.0 

DO  595  JJ  =  1, LONGS 

IF( .NOT. I  NO { JJ)  )  GO  TO  595 

SD=CM$( JJ)-MU 

SS=SS+SO*SD 
595   CONTINUE 

IF(SS.LE.V?S)  RETURN  1 

LA=DSORT(SS) 
LAMBDA  EQUALS  SOUARF  ROOT  PF  SUM  OF  SQUARES 

DO  600  JJ=1, LONGS 

IF(.NOT.INDUJ))  GO  to  *00 
$(I)  NOT  IN  OS  IMPLIES  GAMMA(I)=0 

SD=OM$(  JJ)-MI| 

GAS( JJ) =SD/LA 
GAMMA  VFC.TQR  IS  FTLLED  HERF 
600   CONTINUE 

-RETURN 

END 

SUBROUTINE  AWMCAL 
THIS  P*PT  CALCULATES  THF  H»S  AND  AMM 
IMPLICIT  RF*L*B(  A-H,P-Z,$) 
LOGICAL  SH0PTQtFMLl 
DIMENSION  M(3600)  ,QHfOO) 
COMMON/WORKB/-J(3600),Qf3600) 
COMMON/ XFED1/L,M,LL,LM, ML, MM 
300   DO  305  1=1, LONG 
305   H(I)=0.0 
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341 
342 
350 
360 
37  0 


37  4 
37^ 


380 


DO  370 
JK=J(K) 
IF(Y( JK 

CALL  AR 
DO  360 
DO  350 
I=(LI-1 

TDL=IAP 
IDM=IAR 

lDL$=m 
iDM$=m 

S=DSTNC 

B=DEXP( 
IFCLPN^ 
IFCS.GT 

GO  TO  3 
IFCS.GT 
H(I  >=Hf 
CONTTNU 
CONTIMU 
CONTTNU 
DO  375 
IF(P( T) 
0(1 )=DF 
E(I)=P( 
GO  TO  3 
E(T)=0. 
CONTTNU 
IF( III. 
S=0.0 
DO  380 
S=S+F(T 
AMM  =  S 
RETURN 
END 


K=1,K0 


.EQ.0.01  GO  to  370 

A 

I=LL,LM 

T=ML,MM 

-ISI7F+MT 

(LT-L) 

(MI-M) 

+  1 

+  1 

<IDL$, IDW$! 

S) 

EQ.3< 

3.5) 


-600)  GO  TO  341 
5=0.0 


10.5 
)+R* 


)  R  =  0.0 
Yf  JK) 


=1,L0MG 
CQ.O.O)  GO 
<M-H(I  )) 
)-0(T) 

5 


TIT  374 


X1 

I 

7' 
0 

F 

FO.O)  RFTURN 


T=1,L0NG 

) 


SURRPUTTNF  OMCA?  C 
THIS  POUTINF  CALCMLATFS  THC  OMEGAS 

IMPLICIT  RcAL-«< ft-Htn-Zf$) 

LOGICAL  SHO^totFATL1 

COMMON/ XFFP 1 /L , M t  LL , LM  »  ML  »  MM 
400   IFUJJ.en.l)  GO  TO  420 

DO  410  I=1,L0N^ 

OM(I)=0.0 
410   IM(I)=0 

IF( III.NF.O)  30  TO  415 

DO  490  JJ=1,L0NG 

K=KO 

JK=JJ 

GO  TO  420 
415   DO  490  K=1,K0 

JJ=LCNG 

JK=J(K) 
420   CALL  ARFA 

IM( JK)=1 

DO  480  LI=LL,LM 

DO  470  MT=ML»^M 

I=(LI-1)-ISI7F+MI 

IF(P( I) .PQ.O.0)GO  TO  470 

IDL=IARS(LI-L) 

IDM=IARS(MT-M) 

IDL$=IDL+1 

I0M$=IDM+1 

S=DSTNCE( IDL$,  I0M$) 

B=DEXD(-S) 

I«MLCNG.F0.3600)  GO  TO  441 

IFCS.GT.  •*. 51  3=0.0 

GO  TO  442 

441  IFCS.GT. 10. S)  R=0.0 

442  0M( JK)=CV( JK)+R*E(T ) 
470   CONTINUE 
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48  0   CONTINUE 

IE(  JJJ.EQ.DRETURN 

490   CONTINUE 
RETURN 
END 


SUBROUTINE  AREA 

IMPLICIT  RFAL**( A-H,P-Z,$) 

LOGICAL  SH0DTD,PAIL1 

COMMON/ XFFR1/L,M,LL,LM, ML, MM 

L  =  (  JK-D/ISIZE  +  1 

M=JX-ISIZE=ML-1) 

IF(L.LT.ILO)  30  TO  BIO 


LL=L-IPANGE 

GO  TO  815 

810 

LL=1 

815 

IF(L.GT.IHT) 

LM=L+IPANGE 
GO  TO  *25 

GO 

TO  820 

820 

LM=ISI7F 

8?5 

TF(M.LT.TLO) 

ML=M-I RANGE 
GO  TO  835   . 

30 

TO  830 

830 

ML=1 

835 

IF(M.GT.IHT1 

MM=M+IP ANGF 
GO  TO  R45 

GO 

TO  840 

84  0 

MM=I?I7F 

845 

RETURN 
END 

SUBROUTINE  w| 

?VC 

IMDLICIT  RPA! 

!  £  9 

( A-H,0-Z,S) 

LOGICAL  SH0RTD,1 

FAIL1 

30 

DO  35  K=1,K0 
JK=J(K) 
YYY=YY( JK) 
GAA=GAt JK) 

Y(  JK)=YYY+0*i 

3AA 

IFCY(JK) .GT. 

l.OD- 

YUK)=0.0 

GO  TO  ^5 

34 

IF(Y(JK).LT. 

#qqqqq 0999 90999) 

Y(JK)=1.0 

35 

CONTINUE 

RETURN 

END 

GO  TO  35 


SUBROUTINE  RFEFR 
THIS  ROUTINE  LOADS  OSTNCP  MATRIX 
IMPLICIT  RFAL*8( A-H,0-Z,$) 
LOGICAL  SHORTPtFAILl 
70    DO  75  1=1, B 
00  74  M=ltB 
I$=I-1 
M$=M-1 

-SS=IS*T  $+M$*Mt 
S=DSORT(SS) 

74  DSTNCF(  I,M)=S 

75  CONTINUE 
RETURN 
END 
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SUBROUTINE  CHANGE 
THIS  ROUTINE  CONVERTS  noTiM&L  SOLUTION  FROM  400  CASE  TU 
STARTING  SOLUTION  FOP  B600  C  ?\SE 

IMPLICIT  REAL*S( A-H,0-Z,$) 

INTEGER*^  RTGPOW,BIGCOL 

LOGICAL  SHOOTD,FAIL1 

COMMQN/WORKB/PP(3600) , FILL  (3  600) 

0=0.1 

PREC=10.0**(-5) 

V=PRFC/2**.5 

V2=V*V 

AM=1.0 

LONG=3600 

ISI7.E  =  60 

IRANGE=7 

ILO  =  * 

IHI=53 

ITEP=0 

111=0 

SHORTD=. FALSE. 

FAIL1=. FALSE. 

DO  13  I=l,3*no 

P(I)=1. 0/3600.0 

OM(I)=0.0 

GA(I )=0.0 

IM(I)=0 

J(I)=0 
13    CONTINUE 

00  5010  1=1,400 
FIRST  OETFPy.INE  ROW  AND  COLUMN  OF  INOEX  NUMBER 

TEy,PY=Y(I)/^.0 

NROW=U-1)/20+1 

NCOL=I-20*(NR0W-l) 
NOW  TO  FINP  BASE  SOJAPF  WHICH  IS  TOP  LEFT  SOUARE 

BIGRPW=NROW*3-2 

BTGCCL=NC0L*3-2 

IBASF=PIGC0L+60*(RIGR0W-1) 

YV(IBASE)  =TEMPY 

YYUBASF  +  1)       =tfmpy 

YY(TPASE+?)       =TEMPY 

YY(IBASF+60)     =TEMPY 

YYURASE+61)     =TF^py 

YY(IBASF+6?)    =TEMPY 

YY( IBASE+120)=TEMPY 

YY( IBASE+121 )=TFMPY 

YY( IBASF+1?2)=TFMPY 
5010  CONTINUF 
NOW  TH  FILL  Y  E  P  MATRIX 

DO  5020  1=1,3600 
5020  Y(I )=YY(I) 

RETURN 

END 
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